Demographics and Air Quality by US Counties

Introduction + Purpose

For this project, I wanted to create a linear regression model for different population demographics by US location and respective air quality in that location. This was very interesting to me because I am very interested in environmental discrimination such as environmental racism. Therefore, I wanted to know if factors such as poverty, race, and gender affect the air quality in respective areas because this could potentially point to a symptom of environmental discrimination.

Data Description + Sources

The data for poverty was in the form of percentage of people in poverty by US county. The racial data was organized by white, black, hispanic, asian, and native american percentages by US county. The gender data was male percentages by US county. The air quality data was organized by number of days data was recorded, the number of good, moderate, unhealthy, very unhealthy, and hazardous data.

Because data was only available for all counties for 2010 for all of the variables I wanted in my model, that is the data I chose to use, even though it is not the most recent.

The race and poverty by county data is from the US Census’ Data Mapper tool.

The poverty data is from the US Census’ Small Area Income and Poverty Estimates tool.

The air quality data is from the EPA yearly Air Quality Index Report (https://www.epa.gov/outdoor-air-quality-data/air-quality-index-report).

Data Reformatting/Parsing

The air quality data contains features for the number of days the air quality was recorded and subsequently, the number of good, moderate, unhealthy, and hazardous air quality days. The number of days the quality was recorded for all of the counties is not necessarily the same. Therefore, the features for number of good, moderate, etc. days was converted to a percentage of number of days recorded. This was done in R before database storage.

# load the air quality data into a data frame with unnecessary data
airQuality <- read.csv("/Users/vikaba/Documents/Documents/northeastern/second_year/DS4100/final_project/annual_aqi_by_county_2010.csv")[,c(-3,-7,-11:-19)]
# convert all good, moderate, etc. day counts to percentages of days data was recorded
for (i in 1:nrow(airQuality)) {
  for (j in 4:8) {
    airQuality[i,j] = (airQuality[i,j]/airQuality[i,3])*100
  }
}
airQuality$State <- str_trim(as.character(airQuality$State))
airQuality$County <- str_trim(as.character(airQuality$County))
# export data back with the fixed air quality data
write.csv(airQuality, "AirQualityFixed.csv", row.names = FALSE)

The racial percentages were all in different tables. Before putting the data into the database, the race tables were combined into one table in Excel.

The racial, gender, and poverty tables all had “County” after the county names while the air quality data didn’t. In order to match up the data, “County” was deleted in those values. The poverty data also had state names as numbers and state initials in the county field. The initials were deleted and state numbers were converted to String names. This was all done in R before database storage

# load the race, gender, and poverty data into a data frame
race <- read.csv("/Users/vikaba/Documents/Documents/northeastern/second_year/DS4100/final_project/percent_race.csv")[-1915,] # row 1915 has spanish spelling
gender <- read.csv("/Users/vikaba/Documents/Documents/northeastern/second_year/DS4100/final_project/percent_male.csv")[-1915,] # row 1915 has spanish spelling
poverty <- read.csv("/Users/vikaba/Documents/Documents/northeastern/second_year/DS4100/final_project/percent_poverty.csv")
# delete all rows with states as Puerto Rico because the poverty data does not have Puerto Rico data
gender <- gender[gender$State != "Puerto Rico",]
race <- race[race$State != "Puerto Rico",]
# format race and gender data to not have "County" in the county name field
library(stringr)
gender$State <- as.character(gender$State)
gender$County <- as.character(gender$County)
race$State <- as.character(race$State)
race$County <- as.character(race$County)
for (i in 1:nrow(gender)) {
  gender[i,2] <- gsub("County", "", gender[i,2])
  gender[i,2] <- gsub("Borough", "", gender[i,2])
  gender[i,2] <- str_trim(gsub("Census Area", "", gender[i,2]))
  race[i,2] <- gsub("County", "", race[i,2])
  race[i,2] <- gsub("Borough", "", race[i,2])
  race[i,2] <- str_trim(gsub("Census Area", "", race[i,2]))
}
# format poverty data to not have state initials in county name and "County" and convert state numbers to state name
statesAbbs <- data.frame(state.abb, state.name)
statesAbbs$state.abb <- as.character(state.abb)
statesAbbs$state.name <- as.character(state.name)
poverty$State <- as.character(poverty$State)
poverty$County <- as.character(poverty$County)
for (i in 1:nrow(poverty)) {
  state <- unlist(strsplit(poverty[i,2], "[()]"))[2]
  poverty[i,1] <- statesAbbs[match(state, statesAbbs$state.abb),2]
  poverty[i,2] <- gsub("County", "", poverty[i,2])
  poverty[i,2] <- gsub("Census Area", "", poverty[i,2])
  poverty[i,2] <- gsub("Borough", "", poverty[i,2])
  getRidOfPars <- unlist(strsplit(poverty[i,2], "[(]"))[1]
  poverty[i,2] <- str_trim(getRidOfPars)
}
# export data back with the fixed race, gender, and poverty data
write.csv(race, "RaceFixed.csv", row.names = FALSE)
write.csv(gender, "GenderFixed.csv", row.names = FALSE)
write.csv(poverty, "PovertyFixed.csv", row.names = FALSE)

Database Storage

A relational SQL database was then created in my mySQL to store all values for air quality measurements and the demographic measurements by state.

The database has tables for poverty, race, gender, citizenship, and air quality.

# install and load package to connect to mySQL
#install.packages("RMySQL")
library(RMySQL)
# establish a connection with the local envi_model database
driver <- dbDriver("MySQL")
conn <-
  dbConnect(driver,
            user = "root",
            pass = "1234",
            dbname = "envi_model")

There are multiple values for each state in the air quality data. Therefore, these values were aggregated by averaging in mySQL before joining all of the tables to create a master table.

In order to create a predictive model, the tables in the database are all then joined to create a master table of demographics and air quality by state.

# aggregate the air quality table by state and average all of the values for each of the columns
# join all of existing measuring tables by state
query <-
  dbSendQuery(
    conn,
    statement = "SELECT DISTINCT
    air_quality.state,
    air_quality.county,
    poverty.poverty_percent,
    gender.male_percent,
    race.white_percent,
    race.black_percent,
    race.hispanic_percent,
    race.asian_percent,
    race.native_percent,
    air_quality.good_days,
    air_quality.moderate_days,
    air_quality.unhealthy_days,
    air_quality.very_unhealthy_days,
    air_quality.hazardous_days
FROM
    (SELECT DISTINCT * FROM gender) AS gender,
    (SELECT DISTINCT * FROM race) AS race,
    (SELECT DISTINCT * FROM poverty) AS poverty,
    (SELECT DISTINCT * FROM air_quality) AS air_quality
WHERE
    gender.state = gender.state AND
    gender.state = race.state
        AND gender.county = race.county
        AND race.state = poverty.state
        AND race.county = poverty.county
        AND poverty.state = air_quality.state
        AND poverty.county = air_quality.county
ORDER BY air_quality.state"
  )
allData <- fetch(query, n=-1)
write.csv(allData, "AllData.csv", row.names = FALSE)

Data Visualizations

Before creating a model for and evaluation the data, I used Tableau to create several visualizations that demonstrate some of the features by state and county as the data is interesting to visualize.

As can be seen from this visualization, the smaller circles which represent a low percentage of white people, are mostly darker colored than the bigger circles, meaning they have a higher poverty percentage.

Evaluating the Data

Outliers

First, outliers were identified in the response variable of good days by seeing which values were about 3 standard deviations from the mean.

Because this model deals with air quality data, I did not think it was necessary to omit these outliers. From looking at the outliers (of which there were only 8), the number of good days was low for those days because moderate or unhealthy days were high and I think that is important data to have. It also did not seem like these outliers were due to wrongful experimentation.

goodDaysMean <- mean(allData$good_days) # good days mean
goodDaysStDev <- sd(allData$good_days) # good days st dev
allData$stDevFromMean <- ((allData$good_days - goodDaysMean) / goodDaysStDev)
# find any values in good days that are 3 or more standard deviations from the mean
outliers <- allData[which(abs(allData$stDevFromMean) >= 3),]
outliers

Distribution

The response variable of good days was then analyzed for normal distribution with a histogram.

hist(allData$good_days)

A squared transform was then applied to the good days to normalize the good days variable.

# scatter plot of air quality good days data
hist((allData$good_days)^2)

# transform the average good days by square root
allDataGoodDaysNormalized <- allData[,c(-1,-2,-11:-15)]
allDataGoodDaysNormalized$good_days <- ((allData$good_days)^2)

ANOVA

ANOVA’s were done to see if certain features alone were enough to predict the percentage of good days for a location. These were done for male percentage, white percentage, and poverty percentage.

Then, ANOVA’s were done to see if these features with other features present were a significant predictor variable.

Male Percentage Only

# anova to see if male percentage alone is enough to predict number of good days
aov.genderOnly <- aov(formula = good_days ~ male_percent, data = allDataGoodDaysNormalized)
summary(aov.genderOnly)

               Df    Sum Sq   Mean Sq F value Pr(>F)    
male_percent    1 4.923e+08 492314114   88.19 <2e-16 ***
Residuals    1016 5.672e+09   5582325                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Because the p value is less than 0.05, male percentage alone is enough to predict number of good days.

White Percentage Only

# anova to see if white percentage alone is enough to predict number of good days
aov.whiteOnly <- aov(formula = good_days ~ white_percent, data = allDataGoodDaysNormalized)
summary(aov.whiteOnly)

                Df    Sum Sq   Mean Sq F value   Pr(>F)    
white_percent    1 3.172e+08 317232571   55.13 2.39e-13 ***
Residuals     1016 5.847e+09   5754650                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Because the p value is less than 0.05, white percentage alone is enough to predict number of good days.

Poverty Percentage Only

# anova to see if poverty percentage alone is enough to predict number of good days
aov.povertyOnly <- aov(formula = good_days ~ poverty_percent, data = allDataGoodDaysNormalized)
summary(aov.povertyOnly)

                  Df    Sum Sq  Mean Sq F value Pr(>F)
poverty_percent    1 1.602e+07 16018420   2.647  0.104
Residuals       1016 6.148e+09  6051120               

Because the p value is greater than 0.05, poverty percentage alone is not enough to predict number of good days.

All Percentages to See if Poverty, Male, and White Percentages Significant Predictors

# anova to see if poverty, white, and male percentages with other features are significant predictors
aov.goodDaysAll <- aov(formula = good_days ~ ., data = allDataGoodDaysNormalized)
summary(aov.goodDaysAll)

                   Df    Sum Sq   Mean Sq F value   Pr(>F)    
poverty_percent     1 1.602e+07  16018420   3.094 0.078907 .  
male_percent        1 4.818e+08 481782726  93.043  < 2e-16 ***
white_percent       1 2.070e+08 206963851  39.970 3.87e-10 ***
black_percent       1 1.299e+08 129863906  25.080 6.49e-07 ***
hispanic_percent    1 1.383e+07  13826362   2.670 0.102555    
asian_percent       1 5.903e+07  59025829  11.399 0.000763 ***
native_percent      1 2.665e+07  26653409   5.147 0.023492 *  
Residuals        1010 5.230e+09   5178042                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The p values for male and white percentages are less than 0.05, making them significant predictors. The p value for poverty percentage is greater than 0.05, therefore it is not a significant predictor when all of the other features present

Scatter Plots

Scatter plots were also done to see if there was a linear regression relationship between gender, race, and poverty features and air quality numbers

Male Percentage vs Good Days

scatter.smooth(allDataGoodDaysNormalized$male_percent, allDataGoodDaysNormalized$good_days)

There does not seem to be a strong correlation between good days and percent of males. Good days seem to be constant when compared to percentage of males.

White Percentage vs Good Days

scatter.smooth(allDataGoodDaysNormalized$white_percent, allDataGoodDaysNormalized$good_days)

In this plot, the data seems scattered without a significant correlation between white percentage and good days.

Male Percentage vs Good Days

scatter.smooth(allDataGoodDaysNormalized$poverty_percent, allDataGoodDaysNormalized$good_days)

This plot also shows no significant correlation as the points are very scattered.

Spearman Rank and Pearson Moment Coefficients

Spearman and Pearson coefficients were also used to assess correlation between poverty, race, and gender percentages

Before the tests were performed on these features, an effort was made to normalize them for better correlation coefficient testing. No transforms were made on any of the data. The poverty data was relatively normally distributed. The race and gender data was skewed but did not respond well to transforms.

Pearson Moment Coefficient: percent male and number of good air quality days

cor(allDataGoodDaysNormalized$good_days, allDataGoodDaysNormalized$male_percent, method = "pearson")

[1] 0.2826125

Spearman Rank Coefficient: percent male and number of good air quality days

cor(allDataGoodDaysNormalized$good_days, allDataGoodDaysNormalized$male_percent, method = "spearman")

[1] 0.3291413

Pearson Moment Coefficient: percent white and number of good air quality days

cor(allDataGoodDaysNormalized$good_days, allDataGoodDaysNormalized$white_percent, method = "pearson")

[1] 0.2268606

Spearman Rank Coefficient: percent white and number of good air quality days

cor(allDataGoodDaysNormalized$good_days, allDataGoodDaysNormalized$white_percent, method = "spearman")

[1] 0.24137

Pearson Moment Coefficient: percent poverty and number of good air quality days

cor(allDataGoodDaysNormalized$good_days, allDataGoodDaysNormalized$poverty_percent, method = "pearson")

[1] -0.05097768

Spearman Rank Coefficient: percent poverty and number of good air quality days

cor(allDataGoodDaysNormalized$good_days, allDataGoodDaysNormalized$poverty_percent, method = "spearman")

[1] -0.06818472

For all of the features evaluated, the Spearman Rank and Pearson Moment coefficients are relatively close to each other. This means that outliers do not seem to be having an effect on the data. However, all of the coefficients are very low and significantly below 0.8. Because none of the coefficients are 0, there seems to be some correlation between these features and the number of good air quality days, but not enough to say that there is significant correlation. There is the least correlation between percent poverty and good air quality and the most correlation between percent of males and good air quality. None of these correlations, however, are strong.

The race and gender data is also not normally distributed which may be affecting the correlation coefficients. However, because the Spearman Rank and Pearson Moment coefficients for these 2 features were relatively close, the distribution does not seem to be significantly affecting the correlation coefficients

Building the Model(s)

The desired model for this dataset is a multiple linear regression model with some response variable dealing with the number of days of a certain air quality.

Model with Good Days as Response Variable

The first model has the good days as a response variable.

The data was first randomly split up 50-50 into a training and validation sets.

# training data set: random half of the data
training <- allDataGoodDaysNormalized[sample(nrow(allDataGoodDaysNormalized), nrow(allDataGoodDaysNormalized) / 2),]
# test data set: the other half of the data
test <- allDataGoodDaysNormalized[-c(as.numeric(rownames(training))),]

To create a multiple linear regression model, backward fitting was used and features with p values > 0.05 were removed at each step.

# linear regression model using training data and backward fitting, only keeping p values < 0.05
#goodDaysModel <- lm(formula = good_days ~ ., data = training)
goodDaysModel <- lm(formula = good_days ~ .-(hispanic_percent+poverty_percent+asian_percent+black_percent), data = training)
summary(goodDaysModel)

Call:
lm(formula = good_days ~ . - (hispanic_percent + poverty_percent + 
    asian_percent + black_percent), data = training)

Residuals:
    Min      1Q  Median      3Q     Max 
-6714.4 -1757.0    39.4  1700.9  6561.8 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -11948.969   2911.884  -4.104 4.74e-05 ***
male_percent      281.376     59.463   4.732 2.89e-06 ***
white_percent      41.570      6.594   6.304 6.33e-10 ***
native_percent     70.824     19.852   3.568 0.000395 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2263 on 505 degrees of freedom
Multiple R-squared:  0.1387,    Adjusted R-squared:  0.1336 
F-statistic: 27.12 on 3 and 505 DF,  p-value: 2.823e-16

Next, using the test data set was used to judge the accuracy of the model by calculating the prediction accuracy with the predicted values for the test data set. AIC and BIC calculations were also performed as they are good measures for model selection, in case other models are created for this data.

testUnnormalized <- test
testUnnormalized$good_days <- allData[-c(as.numeric(rownames(training))),c("good_days")]
goodDaysPreds <- predict(goodDaysModel, testUnnormalized)
# data frame of actual values vs predicted values
actuals_preds_good <- cbind(data.frame(actuals = testUnnormalized$good_days, predicteds = goodDaysPreds))
# prediction accuracy
correlation_accuracy <- cor(actuals_preds_good)
paste("Prediction Accuracy: ", correlation_accuracy[1,2] * 100, "%", sep="")

[1] "Prediction Accuracy: 38.6344403240153%"

# fit AIC and BIC calculations
paste("AIC:", AIC(goodDaysModel))

[1] "AIC: 9314.03660913136"

paste("BIC:", BIC(goodDaysModel))

[1] "BIC: 9335.19884921411"

Using the training and test data set to calculate the prediction accuracy, the model was correct about 39% of the time.

Model with Combo of Good and Moderate Days as Response Variable

To create this model, the sum of good and moderate days was taken to create the feature of acceptable air quality days.

allDataSumGoodModerate <- allData[,-15]
# create acceptable days feature
allDataSumGoodModerate$acceptable_days <- allData$good_days + allData$moderate_days
# only leave acceptable days feature in and omit out the rest of the air quality features
allDataSumGoodModerate <- allDataSumGoodModerate[,c(-1,-2,-10:-14)]

Then, outliers were identified again in the same method of standard deviation from the mean as the last model. Once again, I decided to leave these outliers in for more accurate data and because I believe this data is significant for analysis

acceptableDaysMean <- mean(allDataSumGoodModerate$acceptable_days) # acceptable days mean
acceptableDaysStDev <- sd(allDataSumGoodModerate$acceptable_days) # acceptable days st dev
# find any values in good days that are 3 or more standard deviations from the mean
allDataSumGoodModerate$stDevFromMean <- ((allDataSumGoodModerate$acceptable_days - acceptableDaysMean) / acceptableDaysStDev)
allDataSumGoodModerate[which(abs(allDataSumGoodModerate$stDevFromMean) >= 3),]

# omit the standard dev from mean columns for model creation
allDataSumGoodModerate <- allDataSumGoodModerate[,-9]

A histogram was then used to see the distribution of the response variable.

Because none of the transform had a significant effect on the distribution of the heavily right skewed data, the data was left as is.

hist(allDataSumGoodModerate$acceptable_days)

The data was then once again randomly split up 50-50 into training and test data sets.

# training data set: random half of the data
trainingAcceptableDays <- allDataSumGoodModerate[sample(nrow(allDataSumGoodModerate), nrow(allDataSumGoodModerate) / 2),]
# test data set: the other half of the data
testAcceptableDays <- allDataSumGoodModerate[-c(as.numeric(rownames(training))),]

Then, to create a multiple linear regression model, backward fitting was again used and features with p values > 0.05 were removed at each step.

# linear regression model using training data and backward fitting, only keeping p values < 0.05
acceptableDaysModel <- lm(formula = acceptable_days ~ .-(poverty_percent+male_percent), data = trainingAcceptableDays)
summary(acceptableDaysModel)

Call:
lm(formula = acceptable_days ~ . - (poverty_percent + male_percent), 
    data = trainingAcceptableDays)

Residuals:
    Min      1Q  Median      3Q     Max 
-72.472  -0.804   0.877   2.065  17.952 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      45.66604    6.96454   6.557 1.36e-10 ***
white_percent     0.53732    0.07117   7.550 2.06e-13 ***
black_percent     0.54617    0.07372   7.408 5.44e-13 ***
hispanic_percent  0.12156    0.02705   4.493 8.70e-06 ***
asian_percent     0.55904    0.11933   4.685 3.61e-06 ***
native_percent    0.60591    0.08652   7.003 8.04e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.636 on 503 degrees of freedom
Multiple R-squared:  0.1406,    Adjusted R-squared:  0.1321 
F-statistic: 16.46 on 5 and 503 DF,  p-value: 4.642e-15

Next, using the test data set was used to judge the accuracy of the model by calculating the prediction accuracy with the predicted values for the test data set. AIC and BIC calculations were also performed as they are good measures for model selection, in case other models are created for this data.

acceptableDaysPreds <- predict(acceptableDaysModel, testAcceptableDays)
# data frame of actual values vs predicted values
actuals_preds_acceptable <- cbind(data.frame(actuals = testAcceptableDays$acceptable_days, predicteds = acceptableDaysPreds))
# prediction accuracy
correlation_accuracy <- cor(actuals_preds_acceptable)
paste("Prediction Accuracy: ", correlation_accuracy[1,2] * 100, "%", sep="")

[1] "Prediction Accuracy: 24.4240395426979%"

# fit AIC and BIC calculations
paste("AIC:", AIC(acceptableDaysModel))

[1] "AIC: 3212.77970365574"

paste("BIC:", BIC(acceptableDaysModel))

[1] "BIC: 3242.40683977159"

Using the training and test data set to calculate the prediction accuracy, the model was correct about 24% of the time.

Evalulation of Models

Model 1: Good Days as Response Variable

Overall, this model is not a very good measure for predicting the percentage of good air quality days a location will have. Both the multiple R squared and the adjusted R squared are very low and are significantly below 0.7. This means that the variation of good air quality days is not significantly explained by this model. The calculated MAD for this model is 28.66 which is not very close to 0, meaning there is a significant amount of deviation in this model. The model is statistically significant because the p-values for each feature as well as the overall model p-value are well below 0.05.

Model 2: Good+Moderate Days as Response Variable

This model is also not a great measure for predicting the percentage of acceptable days (that is, good and moderate days). The multiple R squared and the adjusted R sqaures are even lower in this model than the previous model. This means that even less variation in the acceptable air quality days is described by this model than the previous model. The calculated MAD for this model, however, is 0.074, which is much closer to 0 than the previous model. Meaning that there is much less deviation in this model than the previous model. This makes sense because outside of the 15 outliers, the percentage of acceptable days was mostly around 90%. This model is also statistically significant with overall and individual feature p-values well below 0.05. The AIC and BIC of this model are also much lower than the previous model, meaning this model is more likely to be the true model and is closer to the truth.

Conclusion

Overall, I hypothesized that building a model with these above demographic features would be a good predictor for air quality is US locations. However, from the two models I have built and the data I aggregated, this does not seem to be the case.

There are a few reasons why this was the result:

• There is not enough air quality data. While there were around 3,000 records from the census demographic data, there were only about 1,000 records for the air quality data. This means that areas that could have contributed to a more successful model did not have air quality recorded.

• The correlation actually does not exist. I do not believe this is the case, because it has been shown that the correlation exists for water quality and the many environmental racism cases in the US.

• The 2010 data is too old and in 2010, air quality was not dependent on any demographics. I don’t, however, think this was the case as 2010 was only 7 years ago.